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Abstract. The spinning of slender viscous jets can be described asymptotically by one-dimen- 
sional models that consist of systems of partial and ordinary differential equations. Whereas the 
well-established string models possess only solutions for certain choices of parameters and set- 
ups, the more sophisticated rod model that can be considered as e-regularized string is generally 
applicable. But containing the slenderness ratio e explicitely in the equations complicates the 
numerical treatment. In this paper we present the first instationary simulations of a rod in a 
rotational spinning process for arbitrary parameter ranges with free and fixed jet end, for which 
the hitherto investigations longed. So we close an existing gap in literature. The numerics is 
based on a finite volume approach with mixed central, up- and down-winded differences, the time 
integration is performed by stiff accurate Radau methods. 
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1. Introduction 

The understanding of jet spinning is of interest in many industrial applications, including for 
example drawing, tapering and spinning of glass and polymer fibers [321 1231 116) and pellet manu- 
facturing 131) . Considering the spinning of highly viscous fluids, the unrestricted motion of an 
instationary jet's center-line is an important feature, as experiments show (see "break-up mode 4" 
by Wong et al. [39]). In the context of slender-body theory there exist two classes of one-dimensional 
models for the numerical simulation of such a jet, string and rod models, [U [Sj [T5j [40] . The string 
models are asymptotic systems of leading order that result from the three-dimensional free bound- 
ary value problems of Newtonian fluid flows in a strict systematic derivation using expansions in the 
slenderness ratio e (e <C 1). They consist of balance laws for mass and linear momentum. The more 
sophisticated rod models also possess an angular momentum balance. The rod models follow from 
the cross-sectional averaging of the underlying three-dimensional balance equations, assuming that 
the displacement field in each cross-section can be expressed in terms of a finite number of vector- 
and tensor-valued quantities. The constitutive elements of a (special) Cosserat rod are a curve 
and a director triad specifying the position (center-line) and the orientation of the cross-sections, 
respectively. The one-dimensional material and geometrical laws that are needed to close the model 
are heuristically motivated. The Cosserat rod model is no asymptotic system of leading order but 
contains the slenderness ratio e explicitely in the angular momentum balance. As the rod reduces 
to a string as e — ¥ 0, the Cosserat rod can be considered as e-regularized string. This regularization 
allows the rod to overcome limitations that the strings have in their applicability, in particular 
when dealing with time-dependencies. In this paper we present the first instationary simulations of 
a Cosserat rod in a rotational spinning process for arbitrary parameter ranges with free and fixed 
jet end, where the string models failed so far. 
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A string model for the jet dynamics was recently deduced in a rigorous slender-body asymptotics 
from the three-dimensional free boundary value problem given by the incompressible Navier-Stokes 
equations, 27 . Accounting for inner viscous transport, surface tension and placing no restrictions 
on either the motion or the shape of the jet's center-line, it generalizes the previously developed 
string models for straight [TU1 EH [H] and curved [T21 [301 I3H] center-lines (for a detailed survey 
of literature see [17]). The numerical results investigating the effects of viscosity, surface tension, 
inertia and gravity on the jet behavior coincide well with the experiments of |39j . However, the 
applicability of the string model turned out to be restricted to certain parameter ranges. Neglecting 
surface tension and gravity, already for jets in a stationary, rotational two-dimensional scenario no 
"physically relevant" solutions exist for ReRb 2 < 1 with Reynolds number Re -1 <C 1 and Rossby 
number Rb <C 1 according to [TTJ [2] . The numerical evidence of this inviscid bound was specified 
analytically in [5]; it is ReRb 2 =3/(2 mii^ \Ai\ 3 ) ~ 1.4 with A, root of the Airy Prime function. The 
restricted applicability / validity results from a non-removable singularity in the model equations due 
to an inconsistency entering with the asymptotically deduced boundary conditions that prescribe 
the jet tangent at the spinning nozzle. This limitation can be overcome by a modification of the 
closure conditions; the boundary condition is omitted in favor of an interface condition that avoids 
the occurrence of the singularity and ensures the regularity of the string quantities. This change 
implies a different string model describing an other jet regime. For gravitational spinning Hlod et al. 
[211 122) distinguish between three compatible disjoint jet regimes, i.e. inertial, viscous-inertial and 
viscous regimes, that they successfully investigated using the string equations with appropriately 
chosen closure conditions. The classification of the regimes is transferable to rotational spinning, 
[5]. But, here the regimes do not cover the whole parameter range. Already for the stationary, 
rotational two-dimensional scenario an existence gap of the string solutions is observed for Re <C 1, 
Rb <C 1, [5]. It is handed over to the instationary simulations that break down for viscous fiber 
jets under very high rotations as they occur in industrial production processes of glass wool [27] . 
When surface tension, aerodynamic forces and temperature-dependent viscosity are included, the 
question of existence and solvability becomes much more difficult or even impossible to answer. In 
non-stationary spinning processes the jet behavior and regime might also change over time. To 
handle this difficulty Hlod [20] investigated a numerical (ad hoc) switching of the closure conditions 
in the simulations. The heuristic approach is motivated by the embedding of the instationary string 
equations into the hyperbolic theory of characteristics under certain assumptions. However, the 
studies remain dissatisfactorily in view of real applications. 

The viscous Cosserat rod theory raises hope to open the parameter ranges of practical interest 
and time-dependencies to simulation and optimization. For the coiling of a viscous jet falling 
onto a rigid substrate Ribe [531 131] proposed a rod model with dynamic center-line that allows 
for stretching, bending and twisting and that is clearly superior to the strings in the application 
of a fluid- mechanical " sewing machine" , see stationary simulations in [36 , 9] and stability analysis 
in [3j2 [28]. Based on these studies and embedded in the special Cosserat theory, we developed 
a modified incompressible rod model for spinning [2] that reduces asymptotically to the string 
equations of [37] for a vanishing slenderness parameter e. It not only covers the string models, 
but also overcomes all thitherto restrictions. In case of stationarity the rod solutions exist for all 
parameter ranges and spinning scenarios without any exceptions, and the existing string solutions 
belonging to the different jet regimes (different closure conditions) are their asymptotic limits as e — > 
0, see convergence results in [5]. Corresponding stationary rod simulations have been successfully 
applied in the study and design of glass wool production processes, [H [26]. The instationary 
rod is described by a system of partial and ordinary differential equations that becomes stiff for 
small e and hence requires a careful numerical treatment. Apart from this structure a further 
numerical challenge lies in the accurate realization of the angular momentum effects which involves 
the conservation of the orthonormal director triad that is attached to the jet's center-line and 
characterizes the orientation of the cross-sections over time. Posing, in favor of a material law for 
the inner forces, a modified Kirchhoff constraint r = eds that relates the jet tangent r and the 
director d3 via the elongation e, the vector-valued angular velocity can be expressed in terms of the 
tangent and the scalar- valued spin (tangential angular speed). So the angular momentum balance 
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becomes scalar-valued and the temporal evolution of the triad redundant, as the other components 
can be computed a posteriori. Motivated from the numerics of elastic Kirchhoff beams Audoly et 
al. [5] just recently developed a discrete geometric Lagrangian approach and performed instationary 
simulations for a jet lay-down (see also [7j). Thereby, they studied the effect of inertia in the angular 
momentum balance. Its neglect simplifies the numerics due to a change of the equations' structure. 
The numerical handling of a free jet end was addressed as open question and topic of future research. 

In this paper, we propose a finite volume approach with mixed central, up- and down-winded 
differences for the instationary viscous rod with free and fixed end in Lagrangian and Eulerian 
parameterization, respectively. The time integration is performed by stiff accurate Radau methods 
taking into account the differential-algebraic character of the system. The rotational tensor asso- 
ciated to the orthonormal director triad is realized using unit quaternions. The approach enables 
the simulation of two-dimensional and three-dimensional rotational spinning for arbitrary (unre- 
stricted) parameter ranges for which the hitherto investigations longed and failed. We deal with 
inflow-outflow set-ups with fixed domain and inflow set-ups with time-dependent enlarging domain 
and discuss the results in comparison to stationary rod [2l [5] and instationary string simulations 
[37J HH1 ISO] , respectively. So this paper closes a gap in existing literature. 

The paper is structured as follows. After a short survey of the special Cosserat theory for 
viscous jets in Section [21 we formulate the instationary rod model in Lagrangian and Eulerian 
parameterizations. For the resulting initial-boundary value problems we develop a finite volume 
approach with Radau time integration in Section [3l In Section [4] we perform numerical simulations 
for the two practically relevant spinning set-ups of enlarging and fixed flow domains and investigate 
the instationary effects. By allowing for the study of all parameter ranges, the rod model shows 
its large potential in view of simulating and optimizing non-stationary three-dimensional rotational 
spinning processes in industrial applications in future. 

2. Special Cosserat theory for viscous jets 

A jet is a slender long body. Because of its geometry with slenderness ratio e (e <C 1), its 
dynamics might be reduced to a one-dimensional description by averaging the underlying balance 
laws over its cross-sections. The procedure is based on the assumption that the displacement field 
in each cross-section can be expressed in terms of a finite number of vector- and tensor-valued 
quantities. The special Cosserat rod theory consists hereby of only two constitutive elements, a 
curve specifying the position and an orthonormal director triad characterizing the orientation of 
the cross-sections, for details see [T|. It represents a general framework that might be applicable 
to all materials and set-ups. The core of the description are physically reasonable one-dimensional 
geometrical and material laws. In this work we use the incompressible viscous rod model derived 
in [2] , whose asymptotic limit as e — > are the string equations of j27j [30] . Since the model equa- 
tions can be formulated in various ways depending on the choice of parameterization/coordinates 
(Lagrangian or Eulerian), basis (invariant, director or outer basis), reference system (fixed or rota- 
tional), set-up (time-dependent or -independent flow domain, acting forces, 2d or 3d), dimensions 
(with dimensions or dimensionless) and so on, we start our introduction with the general invariant 
description of the rod in a Lagrangian parameterization from which all other re-formulations can be 
straightforward computed. In addition, we explicitly state the model formulations that are relevant 
in the considered spinning application and that form the basis for the development of our numerical 
approach, i.e. inflow set-up with enlarging domain (free jet end) in Lagrangian parameterization 
as well as inflow-outflow set-up with fixed domain in Eulerian parameterization. This choice of 
parameterization yields initial-boundary value problems on given computational domains in both 
cases, which facilitates the numerical treatment. 

2.1. General invariant formulation of instationary viscous Cosserat rod model. A special 
Cosserat rod in the three-dimensional Euclidean space E 3 is defined by a curve r : Q — > E 3 and an 
orthonormal director triad {di,d 2 ,d 3 } : Q ->• E 3 with Q = {(a,t) e K 2 | a e [cr a (t),a b (t)}, t > 0}, 
where a addresses a material cross-section (material point) of the rod. The domain of the material 
parameter is chosen to be time-dependent to allow for inflow/outflow boundaries and free end in 
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the Lagrangian description. Considering the dynamics of an incompressible isothermal viscous 
incrtial jet with circular cross-sections and constant mass density, the rod model consists of four 
kinematic and two dynamic equations that are equipped with specific geometrical assumptions 
(shape-preserving incompressibility) and material laws. Its invariant formulation reads [2] 

a t r = v (2.1) 
<9td k = u> x d k 
d t r = d a v 
dtK = d a u> + u> x K 
gA d t v = d(j n + k 

gd t (j a ■ — ) = d a m + r x n + 1 

with 



A 2 

7 (di ®di +d 2 ® d 2 +2d 3 (g) d 3 ), I = -2- 



4-7T 

3fiA -^- • d 3 , m = 3fil ( di ® di + d 2 ® d 2 + -d 3 ® d 3 ) ■ r = ed 3 



and appropriate initial and boundary conditions. Note that system (|2.1[) includes d a r = r and 
c^o-dic = kx dk. The derivatives of the curve r with respect to time and material parameter are the 
velocity v and the tangent field t. Due to the orthonormality of the directors {di,d 2 ,d 3 }, their 
derivatives imply the existence of the angular velocity u> and the generalized curvature k. Assuming 
sufficient regularity, v, r as well as u>, n are related according to the stated compatibility conditions 
(third and fourth equations in (|2.ip ). The dynamic equations are the balance equations for linear 
and angular momentum with external loads k, 1 (body force and body couple line density) coming 
from the considered application. In case of temperature dependencies a corresponding balance can 
be added straightforward, cf. [4]. The curve r is here chosen as the mass-associated center-line. The 
line density qA with constant mass density g as well as the polar moment of inertia I Q refer to the 
referential circular cross-sectional area A Q and are hence time-independent. The incompressibility 
leads to a shrinking of the cross-sections when stretching the body. During the deformation their 
shapes are assumed to be retained. This is incorporated in the geometrical model for the angular 
momentum being linear in oj with dilatation measure e > 0. In consequence the actual cross- 
sectional area and moment of inertia are given by A — A /e and I — I /e 2 '■ The reference area A a 
could be replaced by A in f|2.1[) . which requires the adding of the evolution equation d t {eA) = 0. 
The algebraic relation t = ed 3 that determines the tangent via e represents a modified Kirchhoff 
constraint allowing for extensibility. By its introduction, the normal contact force components 
n • di, n ■ d 2 become Lagrangian multipliers (variables of the system). The tangential contact force 
n • d 3 and the contact couple m are specified by linear material laws in the spatial derivatives of 
the linear and angular velocities (strain rates) with dynamic jet viscosity //, |33[ 134] . Note that for 
the discussion and a better understanding of the geometrical assumptions and material laws, it is 
most convenient to formulate the rod model (|2.ip in the director basis {di,d 2 ,d 3 } as we will do 
later on. Summing up, the variables of the rod model (|2.ip are r, {di,d 2 ,d 3 }, e, k, v, o>, n • di 
and n • d 2 . 

Remark 1 (Impact of modified Kirchhoff constraint). The applied Kirchhoff constraint r = ed 3 
poses a geometric relation between curve and director triad in favor of a material law for the force 
components n • di, n • d 2 . It allows the reduction of the unknown vector-valued angular velocity uj 
to the scalar spin W — uj ■ d 3 , i.e. 

u = -xdt(-) +W-. 

e V e / e 

A respective reformulation of (|2.1[) involves a scalar-valued angular momentum balance and re- 
nounces the evaluation of the director triad, but it changes the clear system structure towards mixed 
derivatives. For this centerline-spin representation Audoly et al. [6] proposed a discrete geometric 
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Lagrangian method (that is inspired by |24] for the analogous centerline- angle representation of an 
elastic Kirchhoff beam), handling the derivatives algorithmically. Thereby, they neglected the influ- 
ence of inertia in the angular momentum balance, cf. Remark^ We stand for system (|2.1[) . since 
its composition of partial and ordinary differential equations is well suited for (standard) finite vol- 
ume schemes and stiff accurate Runge-Kutta methods whose convergence results and performance 
are well-known, see Section^ 

Remark 2 (Model simplifications). In case of negligible inertia gdt(3 ■ ui/e) = and no outer 
couple 1 = in the angular momentum balance, the viscous rod system (|2.1[) with modified Kirchhoff 
constraint reduces to 



d t r = 


V 






d t T = 








gA d t v = 


d a n + k 






d a w = 


e 3 

M 


\d a v ■ ( T x 

e z 




d a M = 


e 5 le 2 


\\rxd a T\\ 2 - 


e 



(r x d a T) 

treating the scalar tangential angular speed / spin W (cf. Remark\T^ and tangential contact couple 
component M as variables. The respective contact force n becomes with its tangential component 
N = SfiAodvV ■ r/e 3 

n = — r <9 CT — ^- d a h ^r((a CT v)^ • d a T)r x d a r 

e e \ e° \ e e° e 

^rx^T-^ ( (dj-^^- ■ d a r)T + ^((^v)^ ■ d a r)(T ■ d a r)T 



Here, z 1 - = z — (z • t)t/||t|| 2 for any arbitrary vector z £ E 3 . 

The string system that is the asymptotic slenderness limit of the rod [S] has the form 

d t r = v 
d t T = d a v 

gA d t v = d a ( ||~jj" T 

with N given as above. 

System (|2.ip is written in a Lagrangian setting. Thereby, the material parameterization might be 
determined up to orientation and a constant by using an arc-length parameterized reference confi- 
guration. Alternatively, any other time-dependent parameterization can be used for the formulation 
of the model, defined via an orientated bijective mapping 

S(;t) : [a a (t),a b (t)] -> [S(a a (t),t), S(a b (t),t)] = [s a {t), s b (t)\, a^S(a,t). 

Assuming sufficient regularity, a scalar convective velocity u and a spatial Jacobian j belong to S: 

d t S(a,t) =u(S(a,t),t), d a S(a,t)=j(a,t)>0, with d s u(S(a,t),t) = —(a,t). 

j 

The re-parameterization of all fields carries convective terms with speed u into (|2.ip . Choosing u = 
implies the material description. Instead of imposing u explicitly, also a constraint can be prescribed 
such that u becomes the associated Lagrangian multiplier and hence an additional unknown of the 
system. The mostly used constraint is the arc-length parameterization of the jet curve for all times, 
yielding an Eulerian setting. Here, e — j coincide due to the Kirchhoff constraint. Moreover, 
d t S(a,t) = u(S(a,t),t) prescribes the rate of change of the arc-length S(a,t) to the material point 
a; e is a measure for the strain and d s u(S(cr,t),t) = (dte/e)(a,t) the corresponding relative strain 
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rate. The Eulerian (spatial) description is certainly the most intuitive one for flow problems and 
allows for the transition to stationary considerations. 

2.2. Rotational spinning process — two relevant set-ups. In rotational spinning processes [26 , 
viscous liquid jets leave small spinning nozzles located on the curved face of a circular cylindrical 
drum rotating about its symmetry axis, cf. Fig. 12.11 At the nozzle, the velocity, cross-sectional area, 
direction and curvature of a jet are prescribed. Starting from an initial length of zero, the extruded 
liquid jet grows and moves due to viscous friction, surface tension and gravity. Also aerodynamic 
forces might act, see e.g. [4]. In this paper we aim at a numerical treatment of the non-stationarity. 
For simplicity we neglect surface tension, aerodynamic forces and temperature dependencies and 
restrict to external loads rising from gravity, i.e. in the Lagrangian setting k = gA ge g and 1 = 
with gravitational acceleration g and direction e g , ||e g || = 1. However, note that once the numerical 
concept is established, the other effects can be easily added as it is already done in the stationary 
considerations of industrial spinning processes in 26, 4. In the following we focus on the numerical 
simulation of two set-ups that are important for the understanding and study of the industrial 
application: 

Set-up A: inflow with enlarging domain (free jet end) 

Set-up B: inflow- outflow with time-independent (fixed) domain 
For the inflow set-up we choose a Lagrangian (material) description, whereas the inflow-outflow set- 
up is formulated in an Eulerian (spatial) setting. Certainly, every set-up could also be formulated 
in the other parameterization, but this yields free boundary value problems. Our choice instead im- 
plies initial-boundary value problems on given computational domains, which makes the numerical 
treatment undeniably easier. 

To describe the spinning process of interest (Fig. 12. ip . we follow [2] and use the reference frame 
that rotates with the drum. Let fl = fte^ with e^ = — e g be the angular frequency of the rotating 
device, then we introduce the rotating outer basis {ai(i), a 2 (t), a 3 (t)} satisfying d t aj = O x aj, 
i = 1,2, 3. This makes the position of the nozzle and the direction of the inflow time-independent, 
but introduces fictitious rotational body forces and couples in the dynamic equations due to inertia. 
We deal with Jl-adapted velocity and angular speed, i.e. V[j = v - (O x r) and = u> — fl. Note 
that we skip the subscript o in the following to facilitate the readability. Moreover, we state the 
model equations in the director basis {di, d 2 , ds} for reasons of the material laws and geometrical 
assumptions, see Notation |3] The rod model for rotational spinning has eight physical parameters: 
jet density g, viscosity /i, length L, diameter D and velocity U at the nozzle as well as drum 
radius R, rotational frequency fi and gravitational acceleration g. These induce five characteristic 
dimensionless numbers: Reynolds number Re = gUR/fj, as ratio between inertia and viscosity, 
Rossby number Rb = U/(£IR) as ratio between inertia and rotation, Froude number Fr = U / y/gR 
as ratio between inertia and gravity as well as I = L/R and e = D/R as length ratios between jet 
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length, nozzle diameter respectively and drum radius. For the subsequent numerical investigations, 
we make the Lagrangian system dimensionless by scaling the quantities with the following reference 
values: 

oo = r = R, v = U, t =v /r , K = l/r , u) = r /v A = ttD 2 /A, k = gA v 2 /r 
n = nA vo/r = n gy^e 2 / (4Re) , m = (iAlvo/(nro) = irgv 2 r^e 4 /(16Re). 

In the Eulerian setting we use alternatively sq = f*o, and consistently uo = vq for the intrinsic 
velocity. 

Notation 3. To an arbitrary vector field z = Y^=i -^ a i = Sj=i z i^i ^ ^ 3 we indicate the co- 
ordinate tuples for the rotating outer basis by z = (zi, 22,23) G K 3 and for the director basis by 
z = (zi, Z2, Z3) G M . The director basis can be transformed into the rotating outer basis by the 
tensor-valued rotation R, i.e. R = ai ® dj = Rij&i ® aj € E 3 ® E 3 wit/i associated orthogonal 
matrix R = (i2y) = (di • aj) e SO^). For the coordinate tuples, z = R • z holds. The cross-product 
a x A e M 3x3 between a vector a e M 3 and a matrix A <E R 3x3 is defined by (a x A) • z = a x (A • z) 
/or aZ/ z 6 R 3 . Moreover, we abbreviate P& = diag(l, 1, fe), k £ R and e; G R 3 , z = 1,2,3 for the 
canonical basis tuples. 

Set-up A: inflow in Lagrangian parameterization. Let Qt — {(cr, t) £ R 2 | (J G (— £(i),0), = 
<, t 6 (0,T]} be the flow domain enlarging over time where Qo = holds initially for t = 0. The 
rod model for the inflow set-up (i.e. inflow at the nozzle a = —£(t) and (stress-)free jet end a = 0) 
reads 

(2.2) 



R • d t r = 


V 




9 t R = 


— w x R 




<9 t ee3 = 


d a v + K X 


v + ee3 x cj 


<9 t ft = 


d a uj + k x 




<9 t v = 


R^ n + 


«xn)+vxaH — ^ R ■ e g 
Fr 


P 2 -^ = 

e 


^ 4 


k x m) + ee 3 x n + Iq 
e^Re 



with Coriolis and centrifugal forces as well as corresponding couples due to the rotating reference 
frame 

2 r, 1 

— R • eo x v ;, 

Rb Rb 2 



ko_ = • e n x v - ^^R • (e Q x (e n x r)) 



l n = P 2 .If w + ^R.e n )x(a, + ^R. % 



cj 1 _ \ 1 d t e^ 
-e X Rb R - en + Rb^ R ' e ° 



e 

and material laws 

1 3 1 

"3 = 3^ (^"3 + KlW 2 - «2«l), m = -^3-P 2 / 3 • {d a u + KX (jj) 

The material laws can be alternatively expressed in terms of the strain rates dte and dtK, since 
dte — d a V3 + K1V2 — K2V1 and dtK = d a u) + k x uj. In particular, they are linear in the strain 
rates. However, to avoid mixed time-space derivatives when plugging the material laws into the 
balance equations we use the stated spatial representation yielding second spatial derivatives in 
the dynamic equations. A strict classification of the whole system (|2.2p is not possible, but it has 
a hyperbolic-parabolic character with ordinary differential equations for curve and rotation group 
(director triad). The boundary conditions are 

r(-£(t),t) = e 2 , R(-£(t),t) = ei ® ei - e 2 <g> e 3 + e 3 <g> e 2 
e(-£(t),t) = 1, K(-l(t), t) = 0, v(-t(t),t) = e 3 , u{-£{t), t) = 
n(0,t)=0, m(0,t) = 0. 
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Set-up B: inflow- outflow in Eulerian parameterization. Let St = {(s,t) G M 2 | s G (0,^), t G 
(0, T], £ > fixed} be the flow domain fixed over time. The rod model for the inflow-outflow 
set-up (i.e. inflow at the nozzle s — and outflow at a prescribed length s = I) reads 



R-d t ? 
d t R 

d s (ue 3 ) 
d t n + d s {un) 
d t A + d s (uA) 

dt(Av) + d s {uAs,) 
P 2 • (d t (A 2 u) + d s {uA 2 u)) 



v — we 3 

— (to — UK) X R 

d s v + «xv + e 3 xw 

d s U) + K X UJ 



— — (3 S n + k x n) + Am x uj ^ — -j^R • e g 



Re 
4 

(9 s m + k x m) 
Re v ; 



16 

e 2 Re 



e 3 x n + If 



with 



Rb 

In = Pa • A 2 



R • en x Am 



Rb 



AR ■ (e n x (e n x r)) 



A 2 w x — R • en 
Rb " 



and material laws 

n,3 = 3Ad s u, 

The boundary conditions for t G [0,T] are 

r(0,t) = e 2 , R(0, t) = ei (8 ei - e 2 ® 
u(0,t) = l, k(0,*) = 0, v(Q,i) = 
n(£,t) = 0, m(£,t) = 0. 

Appropriate initial conditions are specified later on. 



m = -A 2 P 2 / 3 ■ (d s uj + k x uj) 



e 3 + e 3 ® e 2 

= e 3 , u(0,t) = 0, 



(2.3) 



kn 



— A 2 d*uR ■ en 
Rb 



A(0,t) = 1 



The computation of the stated model equations (|2.2[) . (I2.3|) from the general invariant system 
(|2.1[) is straightforward, but lengthy. For more details about the determination we refer to [2] . The 
systems can be easily simplified to 2d. Note that the dimension plays no role for the development 
of the numerical scheme but the reduction to 2d will be used for the simulation of a bench-mark 
test scenario in Section 2] (cf. Fig. 12. II and Eqs. (|A.7|> . (|A.8[) for 2d rotational spinning under neglect 
of gravity). 



Remark 4 (Unit quaternions for rotations). The rotations R G 50(3) can be parameterized, e.g. 
in Euler angles or unit quaternions |25j . We use unit quaternions since this variant offers a very 
elegant way of formulating and computing the evolution equation for R (second equation of (|2.2[) or 
(|2.3p respectively). Define 



R(q) = 



4 



4 



4 



4 



2(<7i<72 ~ go<73) 



2(<?i<3 , 2 + qma) 
2(<?i<2 , 3 - qoq2) 



-4 



4 



4 



4 



with unit quaternions q = (qo,qi,q2,qa), ||q 
respectively, dtR = — (w — UK) x R becomes = ^4(w 



-4(z) ^ 5 



2(<72<73 + go9i) 
|| = 1. then dtR 

- UK) 

Zl 
Z3 





2(9i«3 + Qoq2) 

2(9293 - gogi) 

-4 



4 



4 



4 



o 

-22 
-23 



21 

-33 

2 2 



= — w x R becomes <9tq = -4(u;) ■ q (or 
qj wrf/i skew- symmetric matrix 

z 3 \ 

-2 2 
21 
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3. Numerical scheme 

Finite volume schemes are well-established for the numerical solution of time-dependent partial 
differential equations for various applications |37] . In the following we focus on the inflow problem in 
the Lagrangian parameterization because of the tricky initialization and the handling of the length 
change £(t). The inflow-outflow problem in the Eulerian parameterization where the domain length 
is fixed is comparatively much easier. A respective scheme can be established straightforward in an 
accordant way (cf. [3] for the 2d scenario). To set up the numerical concept for our inflow problem 
we rewrite (|2.2I) in a more convenient formulation, for this purpose we define Ok as the zero vector 
in M fc . We introduce the vector of unknowns 



<p = (rii,n 2 ,e,r, q,K,v,tu) G 



3 if) 



with vj — uj/e. To take account of the differential-algebraic structure of the underlying model, we 
additionally consider the mapping 

z{4>) = (0 2 ,e,r,q,K,v,ra7) € M 19 

that consists of all variables possessing an evolution equation in (|2.2I) . Finite volume schemes are 
based on the integral form of the governing equations that are expressed in terms of flux functions 
and source terms. Therefore, we summarize the constituents with respect to their physical meaning 
and later used numerical approximation. The upper index u, d, c indicates the respective fluxes 
considered for up-, down-winded and central differences: 



3 1 3 1 \ 

v,0 7 ,en7,0 2 , ^^( K i y 2 - k 2 vi), ^^ P i/3 ■ (« x w ) ) 



f d W=(0 13 in 1) ^n 2) 4 



f c (0,a CT h(</>)) 

and 

p(0Ag(0)) 



0l5 'Fl^ 3 'Fl^ Pl / 3 ' 9CT(etn) ) 



Here, h(</>) = (O15, V3, enj) and g{4>) = (O13, V3, ^3, 0, ew) hold. The remaining source terms are 
collected in q(<fi). Due to this dispartment, the system ()2.2[) becomes 

dtz(<f>) = dJ u (<P) + dj d (cb) + d CT f c (0, d a h(<t>)) + P{4>, 9M)) + q(0) (3.4) 

where the closure relations are incorporated. 

Concerning the space discretization we introduce a constant cell size Act and define the number 
N(t) of dynamic cells for the time-dependent jet length tit) by help of the floor function |_-J , 

£(t) 



N(t) 



Act 



<J U+1}/ 2 = - (N(t) - Act, j = 0, . . . , 2N(t), 



where ct;, i — l,...,N(t) denote the cell centers, cf. Fig. 13.21 The cell edge <Jn+i/2 = cr = 
represents the jet end, and [0-1/2,0-3/2] is the closest dynamic cell to the nozzle located at a — —lit) = 
—t. The jet growth is realized by adding new cells at the nozzle, hereby the dynamics of the physical 
quantities is not considered till the cells have completely left the nozzle and become dynamic. Before 
they are treated as static and initialized by the boundary condition at the nozzle. The number of 
static cells depends on the ongoing length increase in the time interval under consideration [t, t+At], 
we use Mit) = Nit + At) — N{t). The introduction of the static cells before/around the nozzle 
allows the adequate initialization of a jet of length tit) < Aa and a stable numerical treatment of 
the temporal evolution. 

The idea is now to integrate (|3.4p over the control volume / cell [o- i _ 1 / 2 , <J i+1 / 2 ], i = 1, . . . , N(t) 
and to set up a differential algebraic system (DAE) in time for the cell averages <f>i of the unknown 



10 



WALTER ARNE, NICOLE MARHEINEKE, ANDREAS MEISTER, AND RAIMUND WEGENER 



£(t) 

^ > N = N(t) 

initialization ai a „ aN 

.\ I ■ I ■ ■ I ■ 



rr N = N(t + At) 

(Jq (J\ (72 CTJV V ' 



<?l/2 CT3/2 &N+1/2 - 

nozzle ^- jet end 

£(t + At) 

Figure 3.2. Spatial discretization of the growing jet with N(t) equally sized 
dynamic cells (centers are marked by black squares). The cell edge a n+x/2 = <t = 
represents the jet end; the complete cell [<Tu2> 03/2] is the closest dynamic one to the 
nozzle at a — —£(t). For the initialization static cells are introduced before/around 
the nozzle (red circles), whose quantities are given by the nozzle conditions. 



quantities 



,!/) : --^Z f ' /! 0M)d<7, i = l,...,N(t). 

'"1-1/2 



Act 



The resulting DAE with Zi(t) = z((f>i(t)) has the form 

d 1 
At Zl ~ Act 



z '' — - — [(fi+1/2 - fji-i/a) + ft+1/2 _ fj— 1/2) + OT+1/2 ~ fi-1/2)] 



+ — p(0, 9lT g(^)da + — / q(0)da. (3.5) 

To express all constituents in terms of the time-dependent <fri(t), we define numerical flux functions 
H u , H d , H c in an up-, down- winded and central manner according to the behavior of the physical 
fluxes, i.e. we apply the upwind-strategy for the convective terms, the downwind-strategy for the 
normal forces and the central approximation for the viscous parts, 

f™ +1/2 «#"(<&, :=P(&) 

f? + i/2«# d (^+i):=f ( Wi) 

<p t + cf) l+1 h(0 i+ i) - h(</>i) 



f^ 1/2 «ff c (&,<fc + i):=f y , ■ n 
i = 1, . . . , N(t) — 1. The integrals that contain the source terms are approximated by means of 

^- £ +1 ' 2 ptf, d„g{<t>)) da « P^i-x^i) := p gW 2 ( ^ ) . i = 2, . . . , tf(t) 

/ 1+1/2 q(</>) da « q(<fc), * = 1, . . . , tf(t). 

AfT •/<r i _ 1/2 



As for the boundaries at nozzle and stress-free jet end, the proposed discretizations make use of the 
respective boundary conditions - collected in tfinoz and 4> en d - in a natural way. We use 

f« r^(U(± \ (d r ^ t df±\ fC fC ( , h(^l) - h((f) nOZ ) \ 

1/2 (0no 2 ), fi/2~ f (01 ), Ti/2 = f l<?W, ^ I at nozzle 

f ]v+i/2 ~ f"(<M, f^+i/2 ~ f d (^end), ^+1/2=° at jet end 
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In accordance we take f° 3/2 p(c/>, d a g{4>)) dcr/Acr « P{(j> n oz,<j>i) f° r the source term in the first control 

1/2 

volume. Inserting the numerical flux functions H = H u + H d + H c and source term discretizations 
into system (|3.5[) . we finally obtain its semi-discrete analogon 

= -^-{H^i, <j> i+1 ) - Hfa-ufc)) + P(&-i, + q(&). (3.6) 
at Act 

The system of DAEs (|3.6[) is of index 2 (according to the definition of [25]). For the time 
integration stiff accurate implicit Runge-Kutta schemes, e.g. Radau Ha methods, are suitable, cf. 
Remark [5] We employ a constant time step At. The resulting nonlinear system of equations is 
solved with a Newton method. 

Remark 5 (Runge-Kutta methods for DAEs). Consider an autonomous differential- algebraic sys- 
tem du/dt — f(u,v), = g(u,v) with initial value (u,v)(to) = (uq.vq). Let a time discretization 
to < t\ < . . . < tM with step size At n = t n +i — t n be given. An appropriate implicit Runge-Kutta 
scheme of level s (with coefficients A = (ay) £l sxs , b=(ij)e R S J has the form 

s 

U n +1 = U n + At n ^ b jf( k j> l j)> V n+1 = v n + (h ~ V n \ . . . \l s - V n )A~ T b, 

3 = 1 

where the levels hi, U are the solutions of the nonlinear system 

n 

kj = u n + At n aijfjkj, lj), Q=g(ki,k), i = l,...,s. 

3 = 1 

A scheme is called stiff accurate if the coefficients satisfy bi — a S i . Radau Ha methods possess this 
property, for s = 1, 2 they are specified by the following Butcher arrays for the coefficients: 



c 


A 


I 


1 






1 



1/3 


5/12 


-1/12 


1 


3/4 


1/4 




3/4 


1/4 



Note that = Y]- as the Runge-Kutta scheme is invariant with respect to autonomization. The 
Radau Ha method for s = 1 is well-known as implicit Euler method. For details on Runge-Kutta 
methods for DAEs see e.g. [18] . 

The space discretization via the finite volume scheme is of first order convergence. We aim at 
an appropriate, stiff accurate time discretization. In the context of differential-algebraic equations, 
the order of convergence p depends on the index and the Runge-Kutta methods often tend to 
loose an order for the variables associated to the algebraic equations. Considering a DAE with 
index 2, the Radau Ha method with s — 2 has p — 3 for the differential variables and p = 2 for 
the algebraic ones, whereas the implicit Euler method has p — 1 for both kind of variables, [18]. 
This theoretical result is confirmed by our numerical tests for the inflow-outflow problem in the 
Eulerian parameterization on the time- independent space interval [0,^] (using an fixed equidistant 
spatial grid with As = t/N where nozzle and outflow are located at the cell edges s = s X j 2 = 
and s = s N+l / 2 = i, respectively), see Figure 1331 For the inflow problem in the Lagrangian 
parameterization, the jet length £(t) and hence the space discretization is time-dependent. There 
is no strict separation of space and time as in an usual semi-discretization. Therefore, it is not 
surprising that the performance of the Radau Ha method with s = 2 differs to the theoretical result. 
We observe a loss of convergence order due to the chosen discretization / initialization at the nozzle 
boundary. In correspondence to the space discretization we obtain here first order convergence in 
time for both Radau variants, Figure I3T41 Consequently, to obtain higher convergence for the inflow 
problem the spatial scheme needs to be modified. 

Remark 6 (Choice of temporal and spatial grid sizes). For the forthcoming numerical simulations 
of the inflow- outflow problem in the Eulerian parameterization the choice of time step and spatial 
grid size follows the CFL-condition with respect to the intrinsic velocity. In the inflow problem in the 
Lagrangian parameterization At and Acr do not need coercively to be coupled, but it turns out that 
they have to be adapted in view of the parameters (Re ; Rb,Fr ) of the problem. Smaller parameters 
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imply in general faster and larger changes in the jet dynamics which require a finer resolution. 
Otherwise it might happen that the used Newton method does not converge. 




dt 



Figure 3.3. Convergence of Radau Ha methods for DAEs of semi-discretized 
inflow-outflow problem in Eulerian parameterization - in consistence with theory. 
Absolute L 2 (0, £)-error for fixed end time T. Left: differential variables. Right: 
algebraic variables (n). Bottom: The intrinsic velocity u that has a special role in 
(|2.3p on first glance behaves like all the other differential variables. 




Figure 3.4. Convergence of Radau Ha methods for DAEs (|3.6p corresponding 
to inflow problem in Lagrangian parameterization with enlarging domain and free 
end. Absolute L 2 (0, £(T))-error for fixed end time T. Left: differential variables. 
Right: algebraic variables. Straight line with slope f indicates first order conver- 
gence. 
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4. Simulation results 




Figure 4.5. Temporal evolution of rod (three depicted times) in comparison to 
steady result of [5] for inflow-outflow set-up in 2d. Top to bottom: r, a, k, u. Jet 
parameters are £ = 1, e = 0.1, Re = 1. Left: Rb = 1 (string is applicable). Right: 
Rb = 0.1 (string fails). Note that for u the scaling of the axes differs. 
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In this section we demonstrate the applicability of the rod model for simulating spinning pro- 
cesses and investigate the relevance of the instationary effects. For this purpose we present numerical 
results for two-dimensional and three-dimensional rotational spinning. The two-dimensional rota- 
tional spinning under the neglect of gravity can be considered as a bench-mark test scenario that 
has been often used in literature - e.g. for the study of strings in [351 El HH [33 120] and of station- 
ary rods in [5]. In this scenario (Fr — > oo) the rotation matrix R can be parameterized in terms 
of a single angle a £ [— n/2,0], see Fig. 12.11 Moreover, the number of variables reduce while the 
initial-boundary value problems keep their characteristic structure. The respective model equations 
for the Set-ups A and B in 2d are stated for completeness in the Appendix, Eqs. (|A.7[) and (IA.8I) . 

We start with the study of the longtime behavior of the instationary rods in comparison to the 
known stationary results of 2, 5 . Therefore, we consider Set-up B, inflow-outflow with fixed domain 
in Eulerian parameterization in 2d. The length ratios are exemplary chosen as i = 1 and e = 0.1. 
For the instationary rod we use the straight jet as initialization at t = 0. The temporal evolution 
of the rod quantities are illustrated by help of three depicted time points in Fig. 14.51 showing two 
different sets of parameters. We clearly observe the convergence of the instationary solutions to 
the stationary ones as time increases, t —¥ oo. The case Re = 1 = Rb (Fig. 14.51 left) lies in the 
parameter regime where also the string model is applicable. As e — > the rod solution coincides 
with the string solution in consistency to the theoretical results (rod-to-string-convergence proof 
in [5]). The other case Re = 1, Rb = 0.1 represents a jet of same viscosity, but exposed to faster 
rotations. In this regime the strings fail. Figure 14.51 shows here instationary jet simulations which 
can be similarly performed for the general three-dimensional inflow-outflow set-up with Fr < oo. 

Time-dependencies play a crucial role for highly viscous jets. Considering the spinning of highly 
viscous jets (Set-up A, inflow with free jet end) the instationary jet center-line is an important 
feature as experiments [39] and corresponding instationary string simulations [TTJ (37J 1SD] show. 
Figure |4~61 illustrates the well-known effect for the rod model. Whereas for inviscid flows (large Re) 
the jet grows along a trajectory that coincides with the stationary jet curve computed for a certain 
length (according to Set-up B) , the center-line for a viscous flow (small and moderate Re) is clearly 
dynamic. However, for long-time it approaches to the stationary jet behavior near the nozzle. So, 
the stationary simulations of the inflow-outflow set-up turn out to be very suitable for studying the 
jet's properties close to the spinning nozzle. This hypothesis was already used for the design of 
production processes of technical textiles in [3J, and it is now confirmed for all parameter regimes. 

We come now to the results for the jet spinning (Set-up A) with arbitrarily chosen parameters 
and free end, whose efficient numerical handling is of main importance and interest for industrial 
applications in future. Figure 14771 presents a growing jet in 2d for two different sets of parameters. 



O t=0.75 
* steady 





O t=0.95 
* steady 



-0.8 



-1 - 



1.2 



1.4 



1.6 



1.2 



1.4 



1.6 



Figure 4.6. Time-dependencies of center-line for viscous jet spinning up to length 
I (Set-up A) in comparison to stationary result computed with same fixed length i 
in Set-up B. Jet parameters in 2d are e = 0.1, Rb = 1. Left: Re = 100 (jet growing 
along stationary curve). Right: Re = 1 (dynamic curve approaching stationary 
behavior at nozzle for t large). 
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Figure 4.7. Growing jet with free end (Set-up A in 2d), e = 0.1. Left: Re = 1, 
Rb = 4 (classical string regime, cf. string results [29]). Right: Re = Rb = 0.1. 
Note that depicted time points are chosen with respect to the dynamics. 



The case Re = 1, Rb = 4 (Fig. 14.71 left) has been already tackled by Panda [29] using the string 
model. We use here the same depicted time points for the rod with e = 0.1, the results are in 
very good agreement as e — > 0. The other case Re = 0.1 = Rb (Fig. 14.71 right) lies outside of 
the applicability regime of the strings. So far, no respective simulations did exist. As expected 




o 


t=2.5 


+ 


t=5 


* 


t=10 




8r 

6 

4- 

2 



-2 

-4- 

-6- 




-2 
r2 



-2 
r2 



o 


t=0.7 


+ 


t=0.9 


* 


1=1 .1 




Figure 4.8. Growing jet with free end, Set-up A in 3d with 2d projections of 
top view (bottom), e = 0.1. Left: Re = 1, Rb = 2, Fr = 2 (classical string regime, 
cf. string results [30]). Right: Re = Rb = Fr = 0.1. Note that depicted time points 
are chosen with respect to the dynamics. 
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the instationary Cosserat rod model opens the full parameter range to the simulation. This also 
holds true for the three-dimensional rotational spinning including gravity. Figure shows the jet 
dynamics with respect to two arbitrary parameter tuples: one is chosen in the (classical) string 
regime, the other one outside of it. The parameter tuple associated to the string regime is taken 
from [30) . the results of rod and string model agree as in 2d. The string model in general fails 
when facing high forces exposed to (highly) viscous jets (small Re, Rb, Fr). The large deformations 
imply boundary layers which cause singularities in the solution and a break-down of theory (model) 
and numerics. The rod model overcomes this limitation since it is an e-regularization of the string. 
This fact is strict in theory for e ^ and holds numerically as long as e is evidently larger than 
zero (machine precision). As for the performance of the simulations, we already mentioned in 
Remark [5] that the applied discretization (At, Act) depends crucially on the considered problem 
parameters. Smaller parameters imply faster, larger changes in the dynamics and higher elongation 
e which require a finer resolution. Theoretically, all parameter settings could be computed by help 
of our proposed scheme, but the simulations are practically restricted to problems with moderate 
elongation e < 50 due to the drastically increasing computational effort. This is no drawback for 
many spinning processes. But for example in industrial processes (like melt-blown) that are driven 
by turbulent air flows much higher elongations are observed. 



The simulation of viscous jet spinning requires the efficient numerical handling of a slender en- 
larging flow domain with free end that is (highly) dynamic due to acting forces. We have proposed 
a finite volume approach with Radau time integration for the viscous Cosserat rod model consisting 
of a system of partial and ordinary differential equations that becomes stiff for small slenderness 
ratio e. The Cosserat rod is an e-regularization of the well-established asymptotic string models 
whose applicability are restricted to certain parameter settings. We have demonstrated the rod's 
superiority and large potential in view of industrial applications by performing instationary simu- 
lations of rotational spinning under gravity for arbitrary parameter ranges, for which the hitherto 
investigations in literature longed and failed. The work has also addressed the open question of a 
free end. 

The established numerical scheme allows the easy incorporation of further practically relevant 
effects, like temperature dependencies and aerodynamic forces. This will be proceeded in future. 
The time and space discretizations depend on the considered problem parameters: small parameters 
cause fast, large deformations and high elongation which require a fine resolution. Due to the 
drastically increasing computational effort the use of our scheme is practically limited to problems 
with moderate elongations, so far. In view of turbulence-driven processes that yield very high 
elongations we intend to get rid of this bottle-neck by investigating appropriate refinement strategies. 



For completeness we state here the model simplifications of (|2.2p and (|2.3[) for the two-dimensional 
rotational spinning scenario where gravity is neglected (Fr — > oo), cf. Fig. 12.11 The rotation is 
parameterized in terms of a single angle a G [— tt/2, 0]. Moreover, we set v = (i> 2 , v 3 ), = (— v 3 , v 2 ), 
u = u>i, and all other quantities analogously. With fl = Cli and rotation matrix 



5. Conclusion 



Appendix 




cos a 



sm a 




we obtain after renumbering, i.e. (22,23) becomes (zi,z<2) for all z, the following two-dimensional 
systems. 
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Set-up A in 2d: inflow in Lagrangian parameterization, cf. 

R(o) • d t r = v (A.7) 
dta = lu 

dt(ee2) — d<j\i + km 1 - + weei 
d t K = d a uj 

d t v = i-(9 CT n + m ±) - ^ - + _L R(a ) . f 



with 



_ u 4 16 1 d t e 



„ 1 ,„ x 3 <9 CT w 

«2 = 3^(d CT W2 + KVl), m= 4 — 3~ 



Set-up B in 2d: inflow- outflow in Eulerian parameterization, cf. 

R(a) ■ d t f = v - we 2 (A.. 

SjOf = CJ — UK 

d s (ue2) = d s v + kv + ujei 
dtK + d s (uK) = d s oj 
d t A + d s {uA) = 

dt(Av) + d s {uAv) = ^-(d s n + kh^) - Auov 1 - - J-Av^ + -L-AR{a) ■ r 
Re Kb Rb z 



d t {A 2 uj) + d s {uA 2 uj) = ^-d s m - 4rr n i + ■rrr^ 2 ^" 
Ke e Ke Kb 

with 



3 

ri2 = 3Ad s u, m = — A 2 d s ui. 

For the simulations in Section 21 the system (|A.8|) is initialized with the straight jet that leaves the 
nozzle perpendicularly, 

n(s,0) = a + l r 2 (s,0)=0 a(s, 0) = 

k(s,0) = u(s,0) = 1 m(s,0) = 

A(s,0) = l Av(s,0) = (0,1) A 2 w(s,0)=0. 
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